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Abstract. The Ginzburg-Landau (GL) equations of superconductivity provide a 
computational model for the study of magnetic flux vortices in type-II superconduc- 
tors. In this article we show through numerical examples and rigorous mathematical 
analysis that the GL model reduces to the frozen-field model when the charge of the 
Cooper pairs (the superconducting charge carriers) goes to zero while the applied 
field stays near the upper critical field. 
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1 Introduction 



Superconducting materials hold great promise for technological applications. Espe- 
cially since the discovery of the so-called high-temperature superconductors in the 
1980s, much research has been devoted to understanding the behavior of these new 
materials. While conventional superconductors require liquid helium (3-4 degrees 
Kelvin) to remain in the superconducting state, high-temperature superconductors 
can be cooled with liquid nitrogen (76 degrees Kelvin) — a clear economic advantage. 
Unfortunately, high-temperature superconductors are ceramic materials, which are 
difficult to manufacture into films and wires, but progress is being made all the time. 

High-temperature superconductors belong to the class of type-II superconduc- 
tors. Unlike type-I superconductors, type-II superconductors can sustain a magnetic 
flux in their interior, but this flux is restricted to quantized amounts — filaments that 
are encircled by a current. The current shields the magnetic flux from the bulk, which 
is perfectly superconducting. The configuration resembles that of a vortex in a fluid, 
and the superconductor is said to be in the vortex state. 
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Figure 1: Phase diagram of a type-II superconductor. 



Figure |l| gives a sketch of the phase diagram of a type-II superconductor in the 
neighborhood of Tc, the critical temperature. The two-dimensional phase space is 
spanned by the temperature T and the (magnitude of the) magnetic field H and is 
roughly divided into three subregions. Each subregion corresponds to a particular 
state: the perfectly superconducting (Meissner) state below the lower critical field 
Hci, where no magnetic field can penetrate the medium; the normal state above the 
upper critical field Hc2, where the superconductor behaves like a normal metal; and 
the intermediate vortex state. Above the critical temperature Tc all superconducting 
properties are lost. 

The vortices, and especially their dynamics, determine the current-carrying ca- 
pabilities of a superconductor. Much effort, both experimental and theoretical, is 
therefore being spent on the study of vortex dynamics and, especially, mechanisms to 
inhibit vortex motion when the superconductor is subject to currents and fields. By 
"pinning" the vortices, one prevents energy dissipation and, hence, loss of supercon- 
ductivity. 

Vortices can be studied computationally at various levels of detail using different 
models. The Ginzburg-Landau (GL) model gives a field (continuum) description that, 
although phenomenological and not based on any microscopic quantum-mechanical 
theory, has been used successfully to study both the dynamics and the structure 
of vortex systems in realistic superconductor configurations If. Figures || and |] 
give two examples of computational results obtained with the GL equations. They 
illustrate both the effectiveness and the difficulties of such calculations. 
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Figure ^ shows a vortex configuration in a two-dimensional cross section of a 
twinned superconducting crystal, which was computed from a steady-state solution 
of the GL equations. The twin boundary (an irregularity in the structure of the 
crystal) is visible in the horizontal line through the center; it acts as a pinning site 
for the vortices. The field is perpendicular to the plane of the cross section, which 
measures 128 x 192 coherence lengths (a characteristic length of the order of microns). 
Each dot corresponds to a vortex intersecting the plane of the cross section; the entire 
configuration has approximately 2,700 vortices. The figure shows the level of detail 
one can achieve with the GL model, given sufficient computing power. At the same 
time, it illustrates the level of computational complexity one faces if one uses the GL 
model. 

Figure ^ shows a series of snapshots of a vortex configuration in three dimensions, 
also computed with the GL model. The objective of this computation was to simulate 
vortex motion through columnar defects and study the potential of the latter as 
pinning sites. The defects are visible as twisted straight lines. The vortices are the 
flexible tube-like structures; they move from one defect to another under the influence 
of external forces. The flgure shows the motion of a vortex that is originally pinned 
on a defect. The vortex develops a loop, the loop peels off, the loop expands in both 
directions in a traveling-wave-like scenario, and gradually the entire vortex transfers 
to the next available defect. 
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Figure 3: Kinking-induced motion of vortices through splayed columnar defects. 



Numerical simulations provide the only way to study vortex dynamics at this 
level of detail. They are an invaluable tool for fundamental research, complementing 
experiment and theory. Numerical simulations of realistic superconductors based on 
the GL model, like the ones illustrated in Figs. and ^ are, however, extremely time 
consuming, and it is desirable to use simpler models whenever possible. Here, we 
focus on the "frozen-field model," which is still a continuum model and the closest 
approximation to the full GL model. In the frozen-field model, the superconducting 
phenomena are decoupled from the electromagnetic field, and the latter is prescribed 
through a vector potential. The frozen-field model is much simpler and has been used 
successfully for numerical simulations of vortex systems 0]. 

In this article, we prove that the frozen-field model is obtained as the asymptotic 
limit of the GL model when the charge of the Cooper pairs (the superconducting 
charge carriers) goes to zero while the applied magnetic field stays near the upper 
critical field. (The upper critical field itself depends on the charge of the Cooper pairs 
and increases as the latter decreases.) Because the temperature is constant in the GL 
model, this limit corresponds to fixing the temperature T and moving up vertically 
through the vortex regime to the curve labeled ifc2 in the phase diagram of Fig. 
The convergence rate is second order in the small parameter. 

For more background on the physics of superconductivity we refer the reader 
to the monograph by Tinkham [Q. The original source for the GL equations of 
superconductivity is . A good introduction to the mathematics of the GL equations 
is 0. The dynamics of the GL equations have been studied by several authors; 
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see 107 i]) m and the references cited therein. The present investigation is closely 
related to the work of Du and Gray 0] . 

Section ^ introduces the Ginzburg-Landau equations, Section ^ contains the 
numerical results and Section ^ the analysis. 



2 The Ginzburg-Landau equations 

In the Ginzburg-Landau theory of superconductivity, the state of a superconducting 
medium is described by a complex scalar- valued order parameter ip and a real vector- 
valued vector potential A. If the state varies with time, a third variable — the electric 
potential — is necessary to fully describe the electromagnetic field. The evolution 
of the state variables is governed by the time-dependent Ginzburg-Landau (TDGL) 
equations, 




0, 



/ I Qj^ \ C C 

a\ V0 - —V xV X A + Js + —V X H = i 

\ c ot J Atc A-k 

where the supercurrent density J s is a nonlinear function of ip and A, 
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These equations are supplemented by the boundary conditions, 

n ■ Js = 0, nx{'VxA) = nxH. 



(2.1) 



(2.2) 



(2.3) 



(2.4) 



Here, H is the applied magnetic field, which we assume to be time independent. The 
constants nis and Qs are the mass and charge, respectively, of a Cooper pair (the 
superconducting charge carriers, also referred to as superelectrons); c is the speed of 
light; and h is Planck's constant divided by 2tt. A Cooper pair is made up of two 
electrons, each with charge — e (e is the elementary charge); hence, Qs is negative, 
Qs = -2e. 

The parameters a and f3 are material parameters; a changes sign at the critical 
temperature Tc, a{T) < for T < Tc (superconducting state) and a{T) > for 
T > Tc (normal state); /? is only weakly temperature dependent and positive for all T. 
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The remaining parameters are cr, the normal state conductivity, and 7, the mobihty 
coefficient. The latter is dimensionless and related to the diffusion coefficient Z), 
7 = h/2msD. 



The boundary conditions (|2.4|) express the fact that superelectrons cannot leave 
the superconductor. Also, if no surface currents are present, the tangential compo- 
nents of the magnetic field must be continuous across the boundary. 

The parameters a and (3 are defined phenomenologically, but they can be ex- 
pressed in terms of measurable quantities, such as the superconducting coherence 
length ^ and the London penetration depth A, 

The coherence length and the London penetration depth define the respective char- 
acteristic length scales for the order parameter and the magnetic induction. Both 
depend on the temperature T and diverge as T approaches the critical temperature 
Tc, because of the factor |a;|~^/^. However, their ratio is, to a good approximation, 
independent of temperature. This ratio is the Ginzburg- Landau parameter, 

K = A/e (2.6) 

In high-Tc superconductors, n is of the order of 50-100. 

The electromagnetic variables are the magnetic induction B, the current den- 
sity J, and the electric field E; they are given in terms of A and by the expressions 

r 1 dA 
B = V X A, J = —V X V X A, E = ^ - V0. (2.7) 

47r c ot 

Equation (|2.2| ) is essentially Ampere's law, J = (c/47r)V x B, where the current J is 
the sum of the supercurrent Jg, the transport current Jt = (c/47r)V x H, and a "nor- 
mal" current J„ = aE (Ohm's law). Hence, the GL model uses a quasistatic version 
of Maxwell's equations, where the time derivative of the electric field is ignored. 



The TDGL equations were first given by Schmid |Tl| in 1966 and subsequently 
derived from the microscopic theory of superconductivity by Gor'kov and Eliash- 



berg |12|. Our notation is the same as in Gor'kov and Kopnin 



The solution of the TDGL equations is not unique. Any solution A, 0) defines 
a family of solutions G^lip, A,(j)) parameterized by a sufficiently smooth function x 
of space and time, 

: A, 0) ^ Le"^'^'/^'^\ A + Vx, - • (2.8) 
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This property is known as gauge invariance; the function x is known as a gauge 
function. Gauge invariance does not affect the physically measurable quantities (the 
magnetic induction B, the magnetization M = B — H, and the current density J). 
Uniqueness requires an additional constraint, which is imposed through a gauge 
choice. The choice of a proper gauge for the TDGL equations has been a subject 
of considerable debate. The choice is a matter of convenience and may depend on the 
problem under investigation. In this article we adopt a gauge in which, at any time, 
the electric potential and the divergence of the vector potential satisfy the identity 

a0+ (c/47r)V- A = (2.9) 

everywhere in the domain, while A is tangential at the boundary. This choice is 
realized by identifying the gauge x with a solution of the linear parabolic equation 

^^_^Ax = fT0 + — V- A, (2.10) 

c dt 4n ^ ^ 4n ' ^ ' 

subject to the condition n ■ Vx = —n ■ A on the boundary. In 0, it was shown that 
the TDGL equations, subject to the constraint (p. 91), define a dynamical system under 
suitable regularity conditions on H. (In the more general case, where H varies not 
only in space but also in time, the TDGL equations define a dynamical process.) This 
dynamical system has a global attractor, which consists of the stationary points of 
the dynamical system and the heteroclinic orbits connecting such stationary points. 
Furthermore, it was shown that every solution on the attractor satisfies the condition 
V ■ A = (and, therefore, also = 0). Thus, in the limit as t — > oo, every solution 
of the TDGL equations satisfies the GL equations in the London gauge. 



2.1 Nondimensional TDGL equations 

In this section, we render the TDGL equations dimensionless by choosing units for 
the independent and dependent variables. Since we are interested in the collective 
behavior of vortices in the bulk of a superconductor in the limit of weak coupling 
{qs — > 0), we take care to choose the units in such a way that they remain of order 
one as — > 0. (We recall that qg is negative, = — 2e.) 

As qs — i> 0, the coherence length ^ remains of order one, while the penetration 
depth A increases like see Eq. ( |2.5| ). This suggests taking the coherence length 

^ as the unit of length. 

To maintain the diffusion coefficient D = h/2'~fms = ^^(T^/lctl)""^ at order one, 
we measure time in units of ''yh/\a\. 
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The real and imaginary parts of the order parameter are conveniently measured 
in units of ipQ = {\a\/PY^'^, which is the value of ip that minimizes the free energy in 
the absence of a field. 

Next, consider the magnetic field. A fundamental quantity in the theory of 
type-II superconductors is the flux quantum $o, 

$o = ^ = 2vr^. (2.11) 

The flux quantum is the unit of magnetic flux carried by a vortex. Together with the 
coherence length and penetration depth, it defines three characteristic field strengths: 
the lower critieal field Hd, the thermodynamieal eritieal field He, and the upper 
critieal field Hc2, 

Below Hci, a superconductor is in the ideal superconducting (Meissner) state, where 
it does not support magnetic flux in the bulk; above Hc2, it is in the normal state, 
where the magnetic flux is distributed uniformly in the bulk; between Hd and Hc2, it 
is in the vortex state, where magnetic flux is quantized in vortex-like configurations 
(see Fig. |lD. The thermodynamieal critical field He is intermediate between Hd and 
Hc2 and is defined by the identity H^/Stc = ^ipQ\a\; H'^/Stt is the energy per unit 
volume associated with a field H, and |'?/'ol'^l minimum condensation energy, 

which is attained when ip = ipo, so these two quantities are in balance when H = He. 

As qs 0, Hd goes to like \qs\, He remains of order one, and Hc2 grows like 
This suggests that we define field strengths in terms of He- In fact, it is 
convenient to absorb a factor so we adopt Hcy/2 or, equivalently, hc/^X\qs\ as 
the unit for the magnetic field strength. 

With the coherence length as the unit of length and Hcy/2 as the unit of field 
strength, it follows that the vector potential is measured in units of ^Hc\/2. Fur- 
thermore, energy densities are measured in units of HI/Att, which is the same as 

Finally, we define the scalar potential in units of (l/7'?/'ofi:|g<i|)(-fff /47r). Notice 
that this unit remains of order one as — >■ 0, because ^Ig^l is of order one. On the 
other hand, the product g^^, which represents an energy density, tends to zero as 
^ 0. (It remains finite on the scale of the penetration depth.) 

Table |l] summarizes the relations between the original variables and their nondi- 
mensional (primed) counterparts. We adopt the latter as the new variables and work 
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until further notice on the nondimensional problem. We omit all primes. 

Table 1: Nondimensionalization. 



Independent variables 


X — ^x' 

t^{lhl\a\)t' 


Dependent variables 


A = {^Hcy/2)A' 

<t>={l/^iPln\q,\){Hl/An)4>' 


Electromagnetic variables 


B = {H,^2)B' 

J = (ci7,V2/47rO J' 

£;=(l/7^2«l'Z.I)(^cV4<)-E' 


Applied field 
Normal conductivity 


H = {H,^2)H' 



The nondimensional TDGL equations are 















OA 






- AA- - 







= 0, 



where 

with the corresponding gauge condition, 

a(/) + V ■ A = 0. 



ij*[V + -A]ij 



(2.13) 
(2.14) 

(2.15) 
(2.16) 



In deriving Eq. ( p.l4| ), we have made use of the gauge condition ( p.l6| ) and the vector 
identity 

AA = -V X V X A + V(V- A). (2.17) 

If Q is the domain occupied by the superconducting material (measured in units of 
^), then Eqs. ( p.l3| )-( 2.16| ) must be satisfied everywhere Q. At the boundary dQ of 
Q, we have the conditions 



n- Js = 0, nx{VxA) = nxH, n ■ A = 0. 
Here, n is the local unit normal vector. 



(2.18) 
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The electromagnetic variables are given by the expressions 

B = VxA, J = VxVxA, E=-dtA-V4>. (2.19) 
The values of the lower and upper critical fields are 

Hci = {2n\iin)-\ Hc2 = n. (2.20) 

2.2 Link variables 

The combination V + (I/k) A plays a fundamental role; we refer to it as the A-gradient 
and write 

Va = V + -A. (2.21) 
The A-gradient defines the A-Laplacian (or "twisted Laplacian"), 

Aa = Va ■ Va = fv + -a) . (2.22) 



K 



The relation between the A-Laplacian and the ordinary Laplacian is most easily 
illustrated by means of the link variables, 

Ux{x,y,z) = expQy" Ax{^,y, z) d^^ , 

Uy{x,y,z) = expQy" Ay{x,ri, z) dr]^ , (2.23) 

Uz{x,y,z) = exp(^^J A,{x,y,() d(^ . 

(We omit the argument t.) The integrals are evaluated with respect to an arbitrary 
reference point. Each Ufj, (/i = x,y,z) is complex valued and unimodular, U* = U~^. 
The vectors A and U may be used interchangeably. With a slight abuse of notation, 
we have 

jj^^{i/K)jA^ VA = f/*Vf/, Aa = U*AU. (2.24) 

3 Numerical solution 



A parallel code for solving Eqs. (p.l3|) - (p.l8|) has been developed as part of a project 
for large-scale simulations of vortex dynamics in superconductors. Details on these 
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simulations and on the code will be published elsewhere; here, we give only a brief 
overview of the numerical methods and the results of numerical simulations showing 
the behavior of the solution as k increases. 

The algorithm uses finite differences on a staggered grid, making all approxi- 
mations accurate to second order in the mesh widths, and an implicit method for 
the time integration, making the algorithm (essentially) unconditionally stable. The 
code, written in C++, has been designed for a multiprocessing environment; it uses 
MPI for message passing. 

We restrict the discussion to rectangular two-dimensional configurations that are 
periodic in one direction and open in the other. The configurations are assumed to 
be infinite in the third, orthogonal direction, which is also the direction of the applied 
magnetic field, H = (0, 0, Hz). 

3.1 Discretization 

Computational grid. The computational grid is uniform, with equal mesh sizes in 
the X and y direction, = hy = h. A vertex on the grid is denoted by Xij = (xj, i/j) 
and is the point of reference for the grid cell shown in Fig. ^. The indices run through 





? € 




A ■ ■ 













Figure 4: Computational grid cell and definition of the discrete variables. 

the values i = 1, . . . ,Nx and j = 1, . . . , A^j^. We assume periodicity in the x direction 
and take the grid so the vertices with j = 1 and j = Ny are located on the open 
boundary of the superconductor. Thus, the size of the domain is = Nx{Ny — l)h'^. 
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Variables. The discrete variables are introduced so that all derivatives are given 
by second-order accurate central-difference approximations. The scalar variables ip 
and (p are defined on the vertices of the grid, 

ipij = ip{xij), (pij = 4>{xij). (3.1) 

(We use the same symbol for the original field and its discrete counterpart.) Vectors 
are defined at the midpoints of the links connecting adjacent vertices, 

Here, and By denote the unit vectors in the x and y direction, respectively. The 
definition of the discrete supercurrent Js is completely analogous. The link variables, 
defined in Eq. (|2.24| ), are obtained from the vector potential, 

[/ ...= Q{i-/K)A^.^ijhx U . . = g(i/K)^t/;ij/*t/_ (3.3) 

They are therefore also defined on the links. Finally, the magnetic induction B, which 
is a vector perpendicular to the plane and given by the curl of the vector potential, 
is defined at the center of a grid cell. 

The definition of the discrete variables is also illustrated in Fig. ^. 

Note that, because of the location of the grid relative to the boundaries, all scalar 
variables, as well as the x components of all vectors {A^, U^, Js,x, and so forth), are 
defined on a iV^. x A^^ grid, whereas the y components of all vectors and the magnetic 
induction are defined on a V^; x {Ny — 1) grid. 



Boundary conditions. We assume periodicity in the x direction, so we need to 
consider the boundary conditions ( |4.7|) only at y = yi and y = y^y- 

The boundary condition for the order parameter, n ■ Va'^P = 0, becomes 

Uy-i,ltlji,2 - = - Uy.i,Ny-lANy-l = 0, (3.5) 

for i = 1, . . . ,Nx. For the vector potential, we require that dyA^ = Hz and Ay is 
constant {Ay = 0) on the boundary. 
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Operators. The gradient of a scalar is a vector and is therefore defined at the 
midpoint of a hnk connecting two adjacent vertices. Thus, 

= + \Ke^) = h-\(f)i+ij - (3.6) 

with an analogous definition for the y component. The gauge-invariant A-gradient 
Va = V + iA is defined in a similar way, with 

(VA^)^;ij = h-\ilJi+ijU:,,ij - ipij). (3.7) 

Thus, the discrete version of the twisted Laplacian Aa is 

(AA^)ij = h'^itpi+ijU^.ij - 2tljij + ilJi_ijU*.i_^j) 

+hy'^{iljij+iUy.ij - + '^ij-iU*.ij^i). (3.8) 

The discrete version of the (normal) Laplacian is defined in the usual way, 

(All)) -J = h~'^{ipi+ij - 2ipij + ipi-ij) + hy'^iipij^i - 2ipij + ipij-i). (3.9) 

The magnetic induction, which is the curl of the vector potential, takes the form 

Bz;i,j = ^x^i^y.i+lj ~ ^y]ij) ~ (A^-iJ + l — Ax-i,j). (3.10) 

We also need the divergence of the vector potential, which is given by 

(V ■ A)^ . = h~\A,,i^, - A,.i_ij) + h-\Ay,i^, - Ay.ij_i). (3.11) 

Algorithm. For numerical purposes, it is useful to treat the TDGL equations ( |4.3| ) 
and ([4.4|) as two separate equations, which are coupled only through certain fields 
and variables. The electromagnetic potentials (p and A are treated as static variables 
in the order parameter equation, 

{dt - {i/K)<f)) ^ - AaV^ - (1 - IVI')^ = 0. (3.12) 

The local nonlinear part of this equation, 

(9i-(i/«:)0)V^-(l-|^|V = O, (3.13) 

is integrated in the simplest possible manner, 

^P,,,it + At) = e-(^/'')<^-^^* {^,,,(t) + At (l - |V',,f ) ^.,,} . (3.14) 
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The nonlocal part, 



dti^ - Aj^^ = 0, 



(3.15) 



is integrated by using a backward Euler method, where the linear equation system is 
solved with the conjugate gradient method. 

The equation for the vector potential, 

adtA- AA- {1/k)Js-V X H = 0, (3.16) 

is linear and depends only indirectly on the order parameter through the supercurrent. 
If we treat the supercurrent as a static variable, we can integrate the equation easily, 
again using the backward Euler method. In the actual implementation, we use the 
fact that the domain is periodic to do a fast Fourier transform in the x direction, 
which leaves us with a tridiagonal system to solve in the y direction. This procedure 
is considerably faster than using an iterative method, such as the conjugate gradient 
method. 



3.2 Numerical results 

We use a rectangular sample, periodic in the x direction and open in the y direction, 
with Nx = Ny = 128. We take = hy = so the sample measures 64 coherence 
lengths in the periodic direction and 63.5 coherence lengths across. (The coherence 
length ^ is defined in Eq. 



First, we considered this system with n = 200 and an applied magnetic field 
Hz = 0.088k. With a relatively large value of k, the surface barrier for vortex entry 
is low, and the system equilibrates relatively fast [l^, |15|. The equilibration required 



5 X 10^ time steps with At = 0.4. The magnetic field produces an almost perfect 
vortex lattice. Figure ^ gives a contour plot of the density of Cooper pairs {ipl"^ at 
equilibrium; the zeros correspond to the centers of the vortices. 

We then started from the configuration of Fig. |^ to find equilibrium configurations 
for other values of k, varying k from Kmin = 40 to Kmax = 800. In this range, the 
ground states are comparable and similar to the one shown in Fig. ^. Since the 
magnetization of a sample is proportional to 1/k^, the vortex density decreases with 
K] below Kmin; the equilibrium state has fewer vortices, and a comparison becomes 
meaningless. Each equilibration required another 3 x 10^ time steps. 
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20 40 60 80 100 120 

Figure 5: Contours of the density of Cooper pairs, IV^P, for a system with n = 200. 

Figure ^ gives the computed values of the quantities 

5^=||V^«-^._||l2, 55. = ^^^^^^^^P^, (3.17) 

\\^z\\l2 

for different values of k. The data show a behavior hke 1/k^ down to k 40. 

Figure ^ shows the average over x of A^^k — ^x.^max a function of y in the bulk 
of the sample, for different values of k. 

The numerical results show that the solution of the TDGL equations converges 
as K increases; in fact, they show quadratic convergence in the small parameter 1/k. 
Given the fact that the Ginzburg-Landau parameter of high-T^ superconducting ma- 
terials is of the order of 50-100, we conclude that the limiting equation is a practical 
alternative in many applications. The question thus becomes: What is the limiting 
equation, and can we confirm the numerical conclusions by rigorous arguments? We 
address this question in the next section. 
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Figure 6: The quantities 5ip (solid squares) and 5Bz (solid discs) for k 
40, 50, 60, 70, 100, 140, 400, 800. The straight lines correspond to behavior. 



4 Asymptotic analysis 



We now return to the TDGL equations ( |2.13| )- (|2.18|) and consider their limit as 
K ^ cxD. These are our standing hypotheses: 



(HI) VL is bounded in R" {n = 2,3), with a sufficiently smooth boundary dQ, for 
example, dQ of class C^'^. 

(H2) The parameters k and a are real and positive. 

(H3) H is independent of time; as a function of position, it satisfies the regularity 
condition H e [^^'^(n)]" for some a G (i, 1). 



(H4) K > 1; 0- = 0(1) and H = 0{k) as k ^ oo. 
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Figure 7: The average < A^^^ - A^,^^,, > vs. y ioi k = 40, 50, 60, 70, 100, 140, 400, 800. 



The assumptions (HI)— (H3) suffice to prove that the TDGL equations define a 
dynamical system in the Hilbert space 

n;i+°'2 = x (4.1) 

see 0. The space W^^°'''^{VL) is continuously imbedded in nL°°(fi), so ip and 

A are bounded and differentiable with square-integrable (generalized) derivatives. 
(H4) is the operative hypothesis for the asymptotic analysis. 

4.1 Mathematical analysis 

Scaling. We begin by scaling the TDGL equations, taking into account the fact 
that we are interested in the limit as ^ {weak coupling), when the applied field 
is near the upper critical field. The scaling is done by means of the dimensionless GL 
parameter k, which grows like 

Since H = 0{k) as k — > oo, we begin by scaling i3" by a factor k, H = kH'. By 
scaling the vector potential by the same factor k, we achieve that the electromagnetic 
variables are all of the same order. 
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The scalar potential is proportional to the charge density of the Cooper pairs, 
which is 0(|gs|) as 0. Hence, k0 remains of order one. This suggests scaling 
by a factor 

Table | summarizes the relation between the current (nondimensional) variables 
and their scaled (primed) counterparts. We adopt the latter as the new variables and 
work until further notice on the scaled problem. We omit all primes. 



Table 2: Scaling. 



Applied Field 


H = kH' 


Dependent variables 


^ = iP' 
A = kA' 

(f) = K-^(j)' 


Electromagnetic variables 


B = kB' 
J — kJ' 
E = kE' 



After scaling, the relevant parameter is k^, rather than k, so we introduce e, 

e = K-^. (4.2) 

The scaled TDGL equations are 

{dt - i£0) - (V + lAfi; - (1 - iV^nV; = 0, (4.3) 

adtA - AA-eJs-V X H = 0, (4.4) 

where 

J, = _l(^*v^/' - ^W*) - I^PA = -53 + iA)^] , (4.5) 

with the corresponding gauge condition, 

£(70 + V • A = 0. (4.6) 

The boundary conditions associated with Eqs. (|4.3| ) and ( [4.4|) are 

n ■ (V + iA)tp = 0, n X (V X A) = n X i?, n ■ A = 0. (4.7) 

The electromagnetic variables are given by the expressions 

B = V X A, J = VxVxA, E = -dtA - eVcj). (4.8) 



Reduction to homogeneous form. Next, we homogenize the problem. Let Aq 
be the (unique) minimizer of the convex quadratic form Qi = 

Q^[A] = J^[{V ■ Af + \V X A- Hf] dx, (4.9) 

on dom((5i) = {A E [W^'^{^1)]"' : n • A = on d^l}. This minimizer satifies the 
boundary-value problem 

VxVxA-Vx// = 0, V-A = infi, (4.10) 

nx {V X A) = nx H, n- A = on dVl, (4.11) 

in the dual of dom(Qi) with respect to the inner product in [L^(fi)]". The map- 
ping H 1-^ Aq is linear, time independent, and continuous from [W^"'^(f2)]" to 
[H/i+".2(Jl)]" The contribution of the vector Aq to the magnetic field is 

So = V X Aq. (4.12) 

We substitute variables, 

A = Ao + eA', (4.13) 

and rewrite the (scaled) TDGL equations ( [4.3D -( ^^ in terms of ip and A' (omitting 
the primes), 

dt^ + ia-\V-{eA))tlj-{V + i{Ao + eA)yij-{l-\^\^)^ = in fi, (4.14) 

adtA - AA - J, = in Q, (4.15) 

where 

Js = -^(^*V^ - ^V^*) - |^|2(Ao + eA), (4.16) 

and 

n ■ V^A = 0, n X (V X A) = 0, n ■ A = on dn. (4.17) 



Functional formulation. We reformulate the system of Eqs. ( [4.14|) -( fi.l71 ) as an 
ordinary differential equation for a vector- valued function u = {ip, A) from the time 
domain (0, oo) to a space of functions on Q, 



The equation is 



u = iij, A) : [0, oo)^ C^ = [L\n)]^ X [L2(^])]-. 



— + Au = fo{u) +6fi{u), 



(4.18) 
(4.19) 
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where A is the hnear operator in associated with the quadratic form Q = Q[ 



u 



on dom((5) = {u 
nonhnear, 



where 



|V^P + a-^ ((V- + |V X A|2)] da;, (4.20) 
{ip, A) G £^ : n ■ A = on dQ}. The functions /o and /i are 
Mu) = {^,{^,A),a-'F,{^,A)), t = 0,l, (4.21) 

2iAo ■ (V^) - |AoP^ + (1 - iV^nV^, (4.22) 
i(l - (T-^)iV ■ A)^ + 2iA ■ (V^A) - (Ao ■ A)^Ij - |A|V, (4.23) 
0, (4.24) 



(Ao + eA). 



(4.25) 



Given any / = {(p,cr ^F) G C^, the equation Au = / is equivalent with the system 
of uncoupled boundary-value problems 



— Atp = if in Q, 
-AA = F in fi, 



n ■ Vif) = on (9n, 

nxA = 0, n-A = on dVl, 



(4.26) 
(4.27) 



in the dual of dom((5) with respect to the inner product in C^. The operator A 
is selfadjoint and positive definite in hence, its fractional powers are well 
defined, they are unbounded if 6' > 0, and dom(74^/^) is a closed linear subspace of 
y^e,2 ^ [|ye,2(^)]2 X [H/e,2^^)]n. Section 1.4]. 



The solution of Eq. ( ^.19| ) depends on we denote it by Ue- We compare Ue with 
the solution uq of the reduced equation 



du 
dt 



+ Au = fo{u). 



(4.28) 



Theorem 4.1 There exists a positive constant C such that 

||u,(t) -no(t)||wiw < C(||n,(0) -uo(0)||wi+",2 +£), t G [0,T]. (4.29) 



Proof. Let i^^j be the ball of radius R centered at the origin in W^"''"'^. Let Ue G Bn 
and uo G -B_r satisfy Eqs. ( [4.19| ) and ( [4.28|) , respectively, with initial data Me(0) and 
Uq{0). The difference v = — uq satisfies the differential equation 



dv 

— +Av = fo{Ue) - /o(%) + efi{Ue) 



(4.30) 
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or, equivalently, the integral equation 

v{t) = e-'MO) + f e~^'~'^^[h{ue) - /o(no) + ds. (4.31) 

Jo 

From the integral equation we obtain the estimate 



X [ll/oK) - k{u,)\\L2+e\\h{uMLA {s) ds. (4.32) 
The operator norms satisfy the inequalities 

||e-*^||wi+-2 < 1, p(i+")/2e-(*'^)^|Ui+-2 < C{t - s)-(i+")/2; (4.33) 

see 



ITj Theorem 1.4.3]. Furthermore, adding and subtracting terms, we have 



Uue)-Uuo) = (2zAo-(V(^,-V'o))-|AoP(^,-^o) 

+ (1 - IV^.P - |^on(^. - V^o) - V^eV^o(^: - Vo). O) ,(4.34) 

where 

||2iAo ■ (V(^^ - Vo))||l2 < 2||Ao||l-||^£ - ^AoIUm 

< C||-?/;e — '?/'o||vyi+«.2 < C'll'We — UQ\\y^i+cc,2, 
|||^Op(V'e -^o)||l2 < C||Ao||ioo||^, - V^oIIl- 

< C\\lp£ — TpoWwi+a,'! < C\\Ue — Mo||wi+".2, 

and the other terms are estimated similarly. Here, C is some (generic) positive con- 
stant, which may depend on H and VL but not on or uq. (In these inequalities we 
have used the continuity of the imbedding of into W^'"^ fl L°°.) The result is 

an inequality of the type 

\\fo{Ue) - /o(^io)||L2 < C\\Ue - Uq\\y^i+^,2, (4.35) 

showing that /q is Lipschitz from W^"'""'^ to C^. 

Using similar estimates, we show that /i is bounded from VV^'^"'^ to so there 
exists a positive constant C such that 

\\h{ue)\W<C. (4.36) 

Combining the estimates ([4.33|) , ([4.35|) , and (|4.36| ) with the inequality ( [4.32|) , we 
conclude that there exist positive constants Ci and C2 such that 

\\v{t)\\y^i+ca < ||t;(0)||i^i+.,2+£Cit(^-")/2^C2 [\t- s)-^'+''^^^ \\vis)\\wi+^,2 ds. (4.37) 

Jo 
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Applying Gronwall's inequality, we obtain the estimate 

||t;(t)||wi+-,2 < C {\\v{0)\\wi+c.,2 + e) , te [o,r], 
for some positive constant C. I 



(4.38) 



It follows from Theorem [4.1| that, if the initial data are such that ||'Ue(0) — 
7x0(0) ||>vi+Q,2 = 0(1) as e I 0, then 



lim = uo 



in C([0,T]; W^+"'2) for any T > 0. In particular, if \\ue{0) - mo(0)|| 



(4.39) 

0{e), 



then the convergence in Eq. ( 4.39 ) is 0{e). 



4.2 Interpretation and final remarks 

It remains to translate the results back in terms of the original variables. We denote 
the solution of the TDGL equations, Eqs. (|2.13| )- (|2.18| ), by ip^, A^, (p^- The variables 
and (p^ are related by the gauge condition acj)^ + V ■ = at all times. Let 



Let Aoo be the solution of the boundary-value problem 

VxVxA-Vxff = 0, V-A = infi, 
nx(VxA) = nxH, n ■ A = on dQ, 



(4.40) 
(4.41) 



and put -Boo = V x A^. The vector A^ and, hence, B^o do not vary with time. Let 
tpoo satisfy the equations 

dtip - AA^tp - (1 - \ip\'^)ip = inn, n- Vao^ = on dn. (4.42) 

Then it follows from Theorem that there exists a positive constant C such that 

B^{t) — Soo||vK«.2 



IIV'K(t) - V^oo(i^)||iyl + ".2 + 

for all t e [0,r], T > 0. 



\H\ 



Pl/a,2 



(n\ I I II-Bk(O) - BooWwc'.i 1 
(0) - V'oo(0)||i^i+<^,2 H h — 



(4.43) 
(4.44) 
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The approximation (^oo, -Boo) is the "frozen-field approximation." Hence, the 
analysis shows that the solution of the TDGL equations converges to the frozen- 
field approximation, uniformly on compact time intervals [0, T] in the topology of 
X as soon as the initial data satisfy the asymptotic estimates 

IIV'k(O) — ^oo(0)||vKi+'».2 = 0(1) and ||Sk(0) — Soo||vk".2 = o{k) as k — >• 00. Under 
slightly sharper conditions we obtain the order of convergence. 



Corollary 4.1 // 



^\ II-Bk(O) — B^\ \wc,2 



^«(0)-^oo(0)||iyi+-2 = O - and " \/^|| " ' =Q - 



as K ^ 00, then 

mt) - ^ooWlk^... + ll^-(^) . ^-11^-- = O (A) ' (4-45) 
uniformly on compact intervals. 

This result explains the numerical results presented in Section ^ 

Remark 1. The asymptotic approximation procedure can be continued to higher 
order, as can be seen from a formal expansion. The equations for the order parameter 
and the vector potential decouple, and at each order one finds first the vector poten- 
tial, then the order parameter. The vector potential satisfies a linear heat equation; 
for example, the first correction beyond Aoo is k~^A, where A satisfies the equation 

- adtA + AA = Q [^:oVa^^oo] . (4.46) 
Remark 2. The analysis given here differs at several points from the analysis of 



Ref. |10[. First, our scaling is slightly different and, we believe, more in tune with 
the physics; second, our regularity assumptions on the applied field are weaker; third, 
our proofs are more direct; and fourth, our results hold in a stronger topology. 
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